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Abstract 

Recent studies of the hydration of micro- and nanoscale solutes have demonstrated a strong 
coupling between hydrophobic, dispersion and electrostatic contributions, a fact not accounted for 
in current implicit solvent models. We present a theoretical formalism which accounts for coupling 
by minimizing the Gibbs free energy with respect to a solvent volume exclusion function. The 
solvent accessible surface is output of our theory. Our method is illustrated with the hydration 
of alkane-assembled solutes on different length scales, and captures the strong sensitivity to the 
particular form of the solute-solvent interactions in agreement with recent computer simulations. 
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Much progress has been made in the last decade in the understanding of hydrophobic sol- 
vation on different length scales 

ELQ 

. Most of this work has been devoted to study solvation 
of purely repulsive, hard sphere-like solutes, while less attention has been given to the influ- 
ence and incorporation of dispersion or electrostatic contributions. Likewise, an entire field 
in the biophysical community has explored electrostatic solvation effects in the absence or 
uncoupled addition to hydrophobic considerations jj. Recently, however, several computer 
simulations have demonstrated a strong coupling between hydrophobicity, solute-solvent dis- 
persion attractions, and electrostatics. For example, a simulation of explicit water between 
paraffin plates revealed that hydrophobie attraction and dewettin, phenomena are strongly 
sensitive to the nature of solute-solvent dispersion interactions 4 . Similarly, simulations 
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of hydrophobic channels [5J, |y] and nanosolutes p| have shown that charged solutes, which 
attract the dipolar solvent due to increasing electric field strength close to the solute surface, 
strongly affect the dewetting behavior and potentials of mean force (pmf). A fully atomistic 
simulation of the folding of the two-domain protein BphC enzyme [8J] further supported 
coupling by showing that the region between the two domains was completely dewetted 
when solvent-solute van der Waals (vdW) and electrostatic interactions were turned off, but 
accommodated 30% of the density of bulk water with the addition of vdW attractions, and 
85-90% with the addition of electrostatics, in accord with experimental results. Finally, Liu 
et al. recently observed a dewetting transition in the collapse of the melittin tetramer which 
was strongly sensitive to the type and location of the hydrophobic residues proving that 
these observations apply to realistic biomolecular systems 

In this letter we propose a continuum description of solvation that explicitly couples 
hydrophobic, dispersion and electrostatic contributions. Similar to the approach of Parker 
et al. in their study of bubble formation at hydrophobic surfaces [h]], we express the Gibbs 
free energy as a functional of the solute cavity shape, the latter given by the volume exclusion 
function of the solvent, and obtain the optimal shape by minimization. This leads to an 
expression similar to the Laplace- Young equation for the geometrical description of capillary 
surfaces h|, but in contrast to existing approaches explicitly includes the inhomogeneous 
distributions of dispersion and electrostatic contributions as well as curvature corrections. 
Geometry-based approaches similar to our formalism exist in related fields, such as the 
Helfrich description of membranes shapes wetting in colloids and granular media |ll 



and electrowetting [12]. We stress that, as opposed to other implicit solvent models 



the solvent accessible surface (SAS) is an output of our theory. This surface encloses the 
optimal solvent accessible volume and should not be confused with the canonical SAS 
which is simply the union of probe-inflated spheres. We begin by verifying that our method 
is able to describe the solvation of small alkanes on molecular scales. We then demonstrate 
that it captures the strong sensitivity of dewetting and hydrophobic hydration to solute- 
solvent interactions on larger scales for a model system of two alkane-assembled spheres. 
In this striking example the strong hydrophobic attraction decreases almost two orders of 
magnitude in units of the thermal energy, ksT, and dewetting is partially or completely 
suppressed when realistic dispersion and electrostatic contributions are included. We expect 
our approach to be particularly useful in solvation studies of proteins where the hydrophobic 
surfaces are highly irregular and laced with hydrophilic units, and a unified description of 
hydration on different length scales is important Q, ^, Q] . 

Let us consider an assembly of solutes with arbitrary shape and composition surrounded 
by a dielectric solvent in a macroscopic volume W. We define a subvolume V empty of solvent 
for which we can assign a volume exclusion function in space given by v(r) =0 for r G V 
and v(r) = 1 elsewhere. We assume that the surface bounding the volume is continuous and 
closed. The absolute volume V and surface S of V can then be expressed as functionals of 
v(f) via V[v] = J w d 3 r [1 — v(f)] and S[v] = j w d 3 r |Vt>(r)|, where V = Vf is the usual 
gradient operator. The density distribution of the solvent is given by p{v{f)) = pov(f), 
where po is the bulk density of the solvent at fixed temperature and pressure. The solutes' 
positions and conformations are fixed. 

We suggest expressing the Gibbs free energy G[v] of the system as a functional of v(r) 
and obtaining the optimal volume and surface via minimization 5G[v]/5v [r\ = 0. We adopt 
the following ansatz for the Gibbs free energy: 

G[v) = PV[v}+ [ d 3 r 7(v)|Vv(f)| + / d 3 r p(v)U(r) 
Jw Jw 

+ I / d 3 r {Vnv)} 2 e(v). (1) 
z Jw 

The first term in (jlj is the energy of creating a cavity in the solvent against the difference 
in bulk pressure between the liquid and vapor, P = Pi — P v . The second term describes 
the energetic cost due to solvent rearrangement close to the cavity surface in terms of a 
coefficient 7. This interfacial energy penalty is thought to be the main driving force for 
hydrophobic phenomena The coefficient 7 is not only a solvent specific quantity but also 



depends on the local topology of the surface Q], i.e., it is a function of the volume exclusion 
function, 7 = j(v(r)). The exact form of this function is not known. For planar macroscopic 
solvent-cavity interfaces 7 is usually identified by the liquid-vapor surface tension, ji Y , of 
the solvent, which we will also employ here. In the following we make a local curvature 
approximation, i.e. we assume that 7 can be expressed solely as a function of the local mean 
curvature of the interface defined by v, ^{vif)) = / y(H(f)), with H(r) = (fti(r) + /t2(r))/2, 
where Ki and k 2 are the two principal curvatures. We then apply the first order curvature 
correction to 7 given by scaled-particle theory fl^ |. the commonly used ansatz to study the 
solvation of hard spheres, arriving at 

1 (H(r))= llv a + ^H(r)), (2) 

where 5 is a constant and positive length expected to be of the order of the solvent particle 

size Q|. The curvature is positive or negative for concave or convex surfaces, respectively. 

Note that this leads to an increased surface tension for concave surfaces, in agreement 
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with the arguments of Nicholls et al. [15] in their study of alkanes. It has been shown by 
simulations that (j2J) predicts the interfacial energy of growing a spherical cavity in water 
rather well for radii > 3A jlfil j. 

The third term in is the total energy of the non-electrostatic solute-solvent interaction 
given a density distribution p{v). The energy U (f) = J2i Uiif) is the sum of the short-ranged 
repulsive and long-ranged (attractive) dispersion interactions between each solute atom 
i and a solvent molecule. Classical solvation studies typically represent U{ as an isotropic 
Lennard- Jones (LJ) potential, U^j(r) = 4e[(<r/r) 12 — (a/r) 6 ], with an energy scale e and a 
length scale a. The importance of treating dispersion interactions independently as opposed 
to absorbing them in to the surface tension term, has been emphasized by Gallicchio et al. 
in their study of cyclic alkanes 

The fourth term in (JTJ) describes the total energy of the electrostatic field expressed by 
the local electrostatic potential ty(v(f)) and the position-dependent dielectric constant e(r) 
assuming linear response of the dielectric solvent. In general, the electrostatic potential 
\I/ can be evaluated by Poisson's equation, V • [e(f)V\l/(f)] = — A(r)/eo, where A(r) is the 
solute's charge density distribution. The most common form for the dielectric function, e(r*), 
is proportional to the volume exclusion function v (r) 

e(v(r}) = e v + v(r)(ei - e v ), (3) 



where e v and e\ are the dielectric constants inside and outside the volume V, respectively. 

Plugging in (J2J) and (jHJ) in functional Q and using the calculus of functional derivatives, 
the minimization yields 

= P - 2 7 i v [H(r) + 5K(r)] - p U(r) - ^[V^(f)e(f)} 2 (---). (4) 

Eq. is an ordinary second order differential equation for the optimal solvent accessible 
volume and surface expressed in terms of pressure, surface curvatures, dispersion interac- 
tions, and electrostatics, all of which have dimensions of force per surface area or energy 
density. K(f) = Ki(r)K2(r) is the Gaussian curvature and follows from the variation of 
the surface integral over H(r) in (|T|). Thus, in our approach the geometry of the surface, 
expressed by H and K, is directly related to the inhomogeneous dispersion and electrostatic 
energy contributions. Note that the SAS is presently defined with respect to the LJ centers 
of the solvent molecules. 

In the following we illustrate solutions of (J1J) in spherical and cylindrical symmetries. 
For a spherical solute (JU) reduces to a function of R, the radius of the solvent accessible 
sphere, H = —1/R and K = 1/R 2 . In cylindrical symmetry the SAS can be expressed by 
a one dimensional shape function r(z), where z is the coordinate on the symmetry axis in 
^-direction and r the radial distance to it. The three-dimensional surface is obtained by 
revolving r(z) around the symmetry axis. We express r = r(t) and z = z(t) as functions 



of the parameter t. The principal curvatures are then given by Ki = —z'j (r\Jr' 2 + z' 2 ) and 
k 2 = (z'r" — z"r')/((r' 2 + z' 2 ) 3 ' 2 ), where the primes indicate the partial derivative with 
respect to t. We solve (JU) and Poisson's equation numerically, using standard forward time 
relaxation schemes. 

We now study the solvation of methane and ethane in water and compare our results 
to the SPC explicit water simulations by Ashbaugh et al. [18], in which the alkanes are 
modeled by neutral LJ spheres. The LJ water-atom parameters are e = 0.8941kJ/mol 
and a = 3.45A for CH 4 , and e = 0.7503kJ/mol and a = 3.47A for CH 3 , and the bond 
length of ethane is 1.53A. We fix the liquid-vapor surface tension for SPC water at 300K 
to 7i v = 65mJ/m 2 [16]. Since we deal with water under ambient conditions the pressure 
term can be neglected and the length 5 remains the only free parameter. For methane 
we can reproduce the simulation solvation energy with a fit 5 = 0.85A. This is in good 
agreement with Huang et al. 16| who measured 5 = 0.76 ± 0.05A for SPC water. Solving 



the cylindrically symmetric problem for the diatomic ethane with the same 5 = 0.85A, we 
obtain a fit-parameter- free G = 11.48kJ/mol, which is only 7% larger than the simulation 
results. Alternatively, the best fit 5 = 0.87A reproduces the simulation energy exactly. This 
is surprisingly good agreement given the crude curvature correction we apply and the fact 
that the large curvature of the system varies locally in space. This supports the validity 
of our continuum approach down to a molecular scale. The curvature and shape functions 
H(z), K(z), and r(z) are plotted in Fig. ^together with the vdW surface and the canonical 
SAS obtained from rolling a probe sphere with a typically chosen radius r p = 1AA over the 
vdW surface [3]. Away from the center of mass \z\ > lA the curvatures follow the expected 
trends H ~ —1/R and K ~ 1/R 2 with R ~ 3.lA for the spherical surfaces. The surface 
resulting from our theory is smaller than the canonical SAS, and is smooth at the center 
of mass [z = 0) where the canonical SAS has a kink. Thus our surface has a smaller mean 
curvature at z = and an almost zero Gaussian curvature, which is typical for a cylinder 
geometry for which one of the principal curvatures is zero. These results may justify the 
use of smooth surfaces in coarse-grained models of closely-packed hydrocarbon surfaces, a 
possibility we will now explore with solvation on larger length scales where dewetting effects 
can occur. 

Let us consider two spherical solutes which we assume to be homogeneously assembled 
of CH 2 groups with a uniform density p=0.024A~ 3 up to a radius R = 15A, defined by 
the maximal distance between a CH2 center and the center of the solute. The CH2-water 
LJ parameters are e = 0.5665kJ/mol and a = 3.536A. Similar ones have been used by 
Huang et al. jj] to study dewetting between paraffin plates. The integration of the CH2- 
water LJ interaction over the volume of a sphere leads yields a 9-3 like potential for the 
interaction between the center of the paraffin sphere and a water molecule [19( . The intrinsic, 
nonelectrostatic solute-solute interaction U ss (s) can be obtained in a similar fashion. The 
solvation of the two solutes is studied for a fixed surface-to-surface distance which we define 
as so = ri2 — 2i? , where r i2 is the solute center-to-center distance. We obtain an effective 
SAS radius of one sphere of about R ~ Rq + 2.4A so that the effective surface-to-surface 
distance is roughly s ~ s — 4.8A. Sine we are also interested in the effects of charging up 
the solutes we place opposite charges ±Ze, where e is the elementary charge, in the center 
or on the edge of the two spheres. 

In the following we focus on a separation distance of sq = 8A to investigate the influence 




FIG. 1: Mean H(z) and Gaussian K(z) curvature and shape function r{z) (solid lines) for ethane. 
The canonical SAS (dashed line) from rolling a probe sphere with radius r p = 1.4A over the vdW 
surface (shaded region) is also shown. 

of different contributions to the energy functional on the shape function, r(z), and the 
curvatures, K[z) and H(z). For sq = 8A, it follows that s ~ 3.2A, such that two water 
molecules could fit between the solutes on the z-axis. We systematically change the solute- 
solute and solute-solvent interactions, as summarized in Tab. I. We begin with only the LJ 
repulsive interactions in system I and then add a curvature correction with 5 = 0.75 A, vdW 
attractions, and sphere-centered charges Z = 4 and Z = 5 in systems II- V, respectively. 
To study the influence of charge location, we shift each charge to the edge of the spheres 
such that they are 8A apart and reduce their magnitude to Z = 1 (system VI). The surface 
tension and dielectric constant of the vapor and liquid are fixed to 7i v = 72m J/m 2 , e v — 1, 
and ei = 78, respectively. 

The results for the curvatures and SAS, defined by r(z), for systems I- VI are shown in 
Fig. |21 Away from the center of mass (|z| > 10A) systems I- VI show very little difference. 



System 


8/k 


vdW attraction 


Z 


W(s )/k B T dewetted 


I 


0.00 


no 





-57.6 


yes 


II 


0.75 


no 





-34.1 


yes 


III 


0.75 


yes 





-6.3 


yes 


IV 


0.75 


yes 


4 


-9.2 


yes 


V 


0.75 


yes 


5 


-5.1 


no 


VI 


0.75 


yes 


1 (oc) 


-1.3 


no 



TABLE I: Studied systems for two alkane-assembled spherical solutes. VF(so) is the inter-solute 
pmf. If r(z = 0) the system is 'dewetted'. In system VI the solutes' charge is located off-center 
(oc) at the solute surface. 



The curvatures are H ~ — 1/R and K ~ 1/R 2 with R ~ 17.4A. Close to the center of mass 
(z ~ 0), however, the influence of changing the parameters is considerable. In system I, 
Eq. (@J reduces to the minimum surface equation H(z) = for z ~ 0. For two adjacent 
spheres the solution of this equation is the catenoid r(z) ~ cosh(z), which features zero 
mean curvature («i and K2 cancel each other) and negative Gaussian curvature. This leads 
to a vapor bubble bridging the solutes. When curvature correction is applied (system II) 
the mean curvature becomes nonzero and positive (concave) at z ~ 0, while the Gaussian 
curvature grows slightly more negative. As a consequence the total enveloping surface area 
becomes larger and the solvent inaccessible volume shrinks, i.e. the value of the shape 
function at z ~ decreases. Turning on solute-solvent dispersion attraction amplifies this 
trend significantly as demonstrated by system III. Mean and Gaussian curvatures increase 
fivefold, showing strongly enhanced concavity, and the volume empty of water decreases 
considerably, expressed by r(z = 0) ~ 10. 7A dropping to r(z = 0) ~ 6.3A. These trends 
continue with the addition of electrostatics in system IV. When the sphere charges are further 
increased from Z = 4 to Z = 5 (system IV-^V), we observe a wetting transition: the bubble 
ruptures and the SAS jumps to the solution for two isolated solutes, where r(z ~ 0) = 0. 
The same holds when going from III to VI, when only one charge, Z — 1, is placed at each of 
the solutes' surfaces. Importantly, this demonstrates that the present formalism captures the 
sensitivity of dewettine: phenomena to specific solvent-solute interactions as demonstrated in 
previous studies □□□□□□ . Note that the SAS at \z\ ~ ±2 A is closer to the solutes in 
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FIG. 2: Mean H (z) and Gaussian K (z) curvatures and shape function r(z) for two alkane-assembled 
solutes of radius Rq = 15 A (shaded region) for systems I- VI. Curvatures are not shown for the 
'wet' systems V and VI. 

VI compared to V due to the proximity of the charge to the interface. Clearly, the observed 
effects, in particular the transition from III to VI, cannot be described by existing solvation 
models, which use the SAS |3[, or effective surface tensions and macroscopic solvent-solute 
contact angles as input. 

The significant change of the SAS with the solute-solvent interaction has a strong impact 
on the pmf, W(sq) = G(so) — G(oo) + U ss (so)- Values of W(so = 8A) are given in Tab. I. 
From system I to VI the total attraction between the solutes decreases almost two orders 
of magnitude. Interestingly, the curvature correction (I— TI) lowers W by a large 23.5/c#T, 
even though R ^> 5. A striking effect occurs when vdW contributions are introduced 
(II — >III) : the inter solute attraction decreases by ~ 28fc#T while the dispersion solute-solute 
potential, U ss (so = 8A), changes by only -QAAksT. Similarly, adding charges of Z = 5 (III 
— > V) at the solutes' centers or Z = 1 (III — > VI) at the solutes' surfaces decreases the 



total attraction by 1.2/cbT and 6k#T, respectively. Note that the total attraction decreases 
although electrostatic attraction has been added between the solutes. The same trends 
have been observed in explicit water simulations of a similar system of charged hydrophobic 
nanosolutes 0. 

These results clearly demonstrate that solvation effects and solvent mediated phenomena 
are not only strongly influenced by solute-solvent interactions, but that these interactions 
are inherently coupled. By including coupling, our formalism captures the balance between 
rydrophobic, dispersive and electrostatic forces which has been observed in previous studies 
DQDQDD but never described in a single theoretical framework. Nonpolar and polar 
coupling is expected to be crucial for a complete characterization of biomolecular solvation. 
The present formalism is only limited by the crude curvature and dielectric descriptions cur- 
rently employed. Future efforts to improve these approximations will be critical to accurately 
describe solvation effects on multiple length scales and for more complicated geometries. 
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